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Abstract 

We performed three-dimensional A^-body/SPH simulations to study how mass resolution and other 
model parameters such as the star formation efficiency parameter, C* and the threshold density affect 
structures of the galactic gaseous/stellar disk in a static galactic potential. We employ 10 6 — 10 7 particles 
to resolve a cold (T < 100 K) and dense (ran > 100 cm -3 ) phase as well as diffuse, hot phases. We found 
that structures of the interstellar medium (ISM) and the distribution of young stars are sensitive to the 
assumed threshold density for star formation, n t h. High-nth models with nth — 100 cm -3 yield clumpy 
multi-phase features in the ISM. Young stars are distributed in a thin disk of which half-mass scale height 
is 10 — 30 pc. In low-nth models with nth = 0.1 cm~ 3 , which is usually employed in cosmological galaxy 
formation simulations, the stellar disk is found to be several times thicker, and the gas disk appears 
smoother than the high-n t h models. A high-resolution simulation with high-nth is necessary to reproduce 
the complex structure of the gas disk. The global properties of the model galaxies in low-nth models, 
such as star formation histories, arc similar to those in the high-n t h models when wc tunc the value of C* 
so that they reproduce the observed relation between surface gas density and surface star formation rate 
density. We however emphasize that high-nth models automatically reproduce the relation, regardless of 
the values of C*. In high-n t h models, the difference in star formation histories is within a factor of two 
for two runs with values of C* which differ by a factor of 15. The ISM structure, phase distribution, and 
distributions of young star forming region are also quite similar between these two. From the analysis of 
the mass flux on phase diagram, we found that the timescale of the flow from the reservoir (jih ~ 1 cm~ 3 ) 
to the star forming regions («h ^ 100 cm -3 ) is about five times as long as the local dynamical time and this 
evolution timescale is independent of the value of C* . The use of a high-n t h criterion for star formation in 
high-resolution simulations makes numerical models fairy insensitive to the modelling of star formation. 
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1. Introduction 

A number of physical processes affect the formation and 
evolution of galaxies. Star formation is among the most 
important processes, not only because it largely deter- 
mines the bulk properties of a galaxy, but also because 
the history of star-formation essentially reflects the for- 
mation history of a galaxy. 

Numerical simulation is a powerful tool to study galaxy 
formation. To compare "simulated" galaxies with ob- 
served ones, it is necessary to follow the dynamics of bary- 
onic matter as well as the assembly of dark matter halos. 
Simulations of galaxy formation are, however, often ham- 
pered by the fact that relevant physics are still poorly 



understood. Numerical resolution is another limiting fac- 
tor. In particular, appropriate physical models of star 
formation should be used in high resolution simulations. 
For example, recent simulations of galaxy formation (e.g., 
Governato et al. 2007), have a spatial resolution of sev- 
eral hundreds pc, with the corresponding mass resolution 
of ~ 10 5 Mq. In such simulations, simple models such 
as an isothermal interstellar medium (ISM) are applied to 
galactic gas disks. In addition, individual giant molecu- 
lar clouds (hereafter GMCs) in galaxies are not resolved 
in current simulations, although GMCs are regarded as 
the site of star formation. Thus one often needs to use 
phenomcnological models, to describe the star formation 
processes, which is called subgrid physics. 
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There are many prescriptions ( "subgrid models" ) of star 
formation used in simulations of galaxy formation with 
coarse resolutions. A commonly used technique is to con- 
vert high-density gas elements to collisionless "star" par- 
ticles (e.g., Katz 1992; Navarro & White 1993; Steinmetz 
& Mueller 1994; Mihos & Hernquist 1994; Katz et al. 
1996; Yepes et al. 1997; Thacker & Couchman 2001; Abadi 
et al. 2003; Kawata & Gibson 2003; Sommer-Larsen 
et al. 2003; Springel & Hernquist 2003; Robertson et al. 
2004; Saitoh & Wada 2004; Okamoto et al. 2005; Stinson 
et al. 2006; Govcrnato et al. 2007; Okamoto et al. 2008). 
Typical criteria to spawn star particles arc as follows (e.g., 
Navarro & White 1993; Katz et al. 1996; Stinson et al. 
2006): (1) the physical density is greater than 0.1 cm~ 3 , 
(2) the temperature is lower than ~ 10000 K, and (3) the 
velocity field is converging. If these three conditions are 
satisfied, 'stars' are then formed at a rate following the 
local Schmidt law. Namely, the local star formation rate 
(SFR), dp*/dt, is assumed to be proportional to the local 
gas density, /3 ga s, and inversely proportional to the local 
dynamical time, td yn ~ 1/ y/ Gp gas : 

dp* _ n Pgas f1 % 

~TT * ~1 ' \ l ) 

OX t-dyn 

where C* is the dimensionless star formation efficiency pa- 
rameter. The value of this parameter is usually calibrated 
by the global star formation properties, the Schmidt- 
Kennicutt relation (Kennicutt 1998; Martin & Kennicutt 
2001). Choosing C* ~ 0.01 reproduces the Schmidt- 
Kennicutt relation in the local universe (e.g., Navarro & 
Steinmetz 2000). However, the threshold density is too 
low which does not correspond to typical densities of the 
neutral hydrogen (HI) and molecular hydrogen (H2) gas 
in real galaxies. 

There are several models that assumes higher density 
regions as star forming regions. For instance, Kravtsov 
(2003) adopted ran > 50 cm -3 as the star forming re- 
gions in a cosmological simulation of galaxy formation. 
When he chose a plausible star formation time, t s { : 
dp*/dt = PgasAsf, which is set to constant (= 4 Gyr) in 
his simulation, the Schmidt-Kennicutt relation is also re- 
produced. More recently, Taskcr & Bryan (2006; 2007) 
performed adaptive mesh refinement simulations of the 
ISM in a static halo potential. Their simulations resolve 
individual star forming regions (the minimum cell sizes are 
25 — 50 pc). They compared two variants of star forma- 
tion criteria: (a) ra H > 10 3 cm~ 3 , T < 10 3 K, and C* = 0.5, 
and (b) n H > 0.02 cm- 3 , T < 10 4 K, and C* = 0.05. Both 
models also employ converging flows as one of the star for- 
mation criteria and the star formation law described by 
equation (1). Interestingly, it is found that both the mod- 
els reproduce the Schmidt-Kennicutt relation. The star 
formation histories arc also found to be similar. Therefore, 
their results appear to imply that the global star forma- 
tion properties are not sensitive to the details of star for- 
mation prescriptions. 

In this paper, we examine how numerical prescriptions 
of star formation affect structure of the ISM, and the 
global star formation history (SFH) of a galactic disk. 



In particular, we focus on the threshold density in star 
formation criteria (nth) and the star formation efficiency 
(C*). We adopt two values of density threshold: 0.1 cm -3 
(low-nth model) and 100 cm -3 (high-nth model). We also 
test the effect of star formation efficiency parameter C* in 
high-nth models. We perform high-resolution SPH simula- 
tions (number of SPH particles are 10 6 — 10 7 ) in a galactic 
potential and we compare structure of the ISM and stellar 
disks. We show that, while both models can exhibit simi- 
lar SFHs and the relation of surface gas density to surface 
SFR, only high-rath models have the complex, inhomogc- 
neous, and multiphase ISM. The ISM has a log-normal like 
probability density distribution (PDF) , while a model that 
does not include either star formation and supernova (SN) 
feedback has a power-law like PDF; Star formation and SN 
feedback distort PDF. Interestingly, structure of the ISM, 
stellar disks, and SFRs (SFHs) arc not simply propor- 
tional to C* in high-rath models. This is because the mass 
supply timescale from from the reservoir (ran ~ 1 cm -3 ) 
to the star forming regions (ran ^100 cm ) is ~ 5 £dyn(nR~) 
and this timescale is not affected by the adopted value of 
C* . 

The plan of this paper is as follows. We mention the re- 
lation between numerical resolution and expressible phase 
of the ISM in Section 2. In Section 3, we describe prop- 
erties of a model galaxy and our numerical methods. Our 
results are presented in Section 4. Summary and discus- 
sions are given in Section 5. 

2. Estimate of required resolutions to appropri- 
ately model the star formation 

We consider regions with the density greater than 
100 cm -3 as "star forming regions" . From the Jeans con- 
dition, we can estimate the resolution necessary to ex- 
press the gravitational collapse for a fluid with a given 
density and temperature. In smoothed particle hydrody- 
namics (SPH) simulations, this condition is expressed as 
Mjcans > Anb x ^sph, where Mje ana is the Jeans mass, 
A„b is a typical number of neighbor particles, and raispH 
is the mass of an SPH particle (Bate & Burkcrt 1997; Bate 
et al. 2003; Hubber et al. 2006). The Jeans condition then 
determines resolved regions in phase (p — T) plane. 

Figure la shows the limits for several different particle 
masses. We can accurately treat the gravitational frag- 
mentation of fluid in the upper region of each line in the 
phase diagram . The line is given in equation (3) in Saitoh 
et al. (2006). The curve indicates the thermal balance be- 
tween the cooling and heating that are both adopted in 
this paper (see §3). We assume A n b = 32. It is clear that 
a very higher resolution is required in order to resolve the 
gravitational collapse to dense clouds, e.g., GMCs. 

Figure lb illustrates Mj eans and tosph for a given crit- 
ical density, ra cr i t , where ra cr ;t is defined as values of the 
density ra cr i t at intersection between equilibrium p—T re- 
lation in figure la and constant Mj eans line. The solid 
curve indicates Mj eans as a function of n cr ;t, while the 
dashed curve represents mspH as a function of ra cr ; t . The 
red region overplotted on the dashed curve represents n cr it 
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of simulations in this paper, where the mass range of SPH 
particles in simulations is on the order of 10 2-3 M© (see 
Table 2). For comparison, we choose several simulations 
of galaxy formation (tosph ~ 10 6 Mq in Abadi et al. 2003; 
"t-Sph ~ 10 5 Mq in Governato et al. 2007) and plot the 
corresponding range of the critical density as the blue re- 
gion. 

3. Methods and models 

We investigate 3-D evolution of a gas disk in a static 
disk-halo potential. We assume a Navarro-Frenk- White 
(NFW) density profile (Navarro et al. 1997) for a dark 
matter halo, and a Miyamoto-Nagai model (Miyamoto & 
Nagai 1975) for a stellar disk. For the halo model, we 
adopt a cosmological model of a standard A CDM universe 
(Spcrgcl et al. 2003). The cosmological parameters are 
flu = 0.3, Q A = 0.7, and H = 70 km s" 1 Mpc" 1 . We 
use these parameters to model the dark matter halo. By 
adopting the static potentials for a dark matter halo and 
a stellar disk, we can prevent global disk instabilities and 
artificial disk heating due to the scattering of gas and 
star particles by massive dark matter particles, which are 
inevitable if we use a live halo with low mass resolution. 
We will study the effect of live halo and also the effect 
of the mass resolution of halo particles in forthcoming 
papers. 

Gravitational forces are computed by a parallel tree 
code ASURA (Saitoh et al., in preparation) that uti- 
lizes the special purpose hardware GRAPE. The par- 
allel implementation with GRAPE is based on that of 
Makino (2004). We use an opening angle 9 = 0.5 for a 
cell opening criterion. We only use a monopolc moment. 
Hydrodynamics is followed by the standard SPH method 
(e.g., Lucy 1977; Gingold & Monaghan 1977; Monaghan & 
Lattanzio 1985; Monaghan 1992). The kernel size of each 
SPH particle is determined by imposing the number of 
neighbors to be 32 ±2. We use a cooling function for a gas 
with the solar mctallicity for a temperature from 10 K to 
10 8 K (Spaans & Norman 1997). An uniform heating from 
the far-ultraviolet (FUV) radiation observed in the solar 
neighborhood (Wolfire et al. 1995) is included. We do not 
include heating from an ultra-violet (UV) background ra- 
diation. This is because we consider the structure of the 
ISM in the current environment of the Milkyway galaxy 
in this paper. When we consider the detailed structure of 
the ISM in the galaxy formation process, a careful treat- 
ment of the UV background and the local FUV radiations 
is required since the fluxes are quite stronger than those 
in the local universe in high redshift universe. This is 
beyond the scope of this paper. 

3.1. Halo+disk galaxy model 

We assume that the dark matter density profile is de- 
scribed by a NFW profile: 

PhBlo(x) = ^lf^2 > x = r / r s, (2) 
cnfw = r v i T /r B , (3) 
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Fig. 1. Required mass resolutions to resolve cold and dense 
gas phase (a) and the critical density as a function of Jeans 
mass in fixed cooling and heating rates using in our simula- 
tions (b). The red lines in panel (a) indicate the Jeans limits 
(the definition is shown in the text) with mgpn = 10 6 , 10 3 , 10 2 , 
and 10 Mq (from top to bottom), where mgpg is the mass 
of an SPH particle and we assume that the Jeans mass is 
32 X mgpH- The black curve represents the equilibrium tem- 
perature of the ISM for the adopted cooling and heating with 
the solar abundance, T eq . For details of the cooling and heat- 
ing rates, see in §3.1. Panel (b) shows the critical density, 
"crit i where n cr i t is defined as values of the density ra cr it at 
intersection between equilibrium p — T relation in panel (a) 
and constant Mj eans line, as a function of mass. Thick solid 
curve shows n cr - lt — Mj eans relation and thin dashed curve 
represents n cr ; t — mgpjj relation. We assume the mass of 
each SPH particle is the same. The red region overplot- 
ted on the dashed curve represents n cr ; t of simulations with 
mgpH = 10 2— 3 Mq (this study), while the blue regions over- 
plotted on the solid curve represents mass ranges of several 
simulations of galaxy formation (rngpn ~ 10 B Mq in Abadi 
et al. 2003; m S PH ~ 10 5 Mq in Governato et al. 2007). 
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where p c is the characteristic density of the profile, r is a 
distance from the center of the halo, r s is a scale radius 
for the profile, p cr is the critical density of the universe, 
and A v j r is the virial overdensity (we employ A vir = 340) . 
The halo mass is set to be M vir = 10 12 Mq and the con- 
centration parameter is set to be cnfw = 12 (Klypin ct al. 
2002). Then the profile has r vir = 258 kpc, r B = 21.5 kpc, 
and p c = 4.87 x fO 6 M Q kpc" 3 . 

The stellar disk is assumed to follow the Miyamoto- 
Nagai model: 



model". The other is n th = 0.1 cm" 3 and T th = 15000 K. 
We dub this model the "low-nth model". The low-nth 
model is similar to what is used in previous simulations 
of galaxy formation (e.g., Navarro & White 1993; Katz 
et al. 1996; Thacker & Couchman 2001; Okamoto et al. 
2005; Governato et al. 2007; Okamoto et al. 2008). 

When an i-th gas particle is eligible to form stars, we 
compute the probability psF,i of the particle to spawn a 
new star particle with mass m* iSpa wn during a time-step 
width dt as 

( r dt 
expl — C* 



p*(R,z) 



Mai 

4-7T 



(5) 



where 



R*R 2 + (R* + 3V^ + z* 2 ) (R* + Vz 2 + z* 



[Rl + (R* + sfz^- 



^7) 2 ] 5/2 (z 2 + z* 



.3/2 



with mass M», radial scale length i?„ and vertical scale 
length z*, respectively. R and z are the cylindrical galact- 
centric radius and the height, respectively. Numerical val- 
ues of model parameters are given in Table 1. 

Our model of the gaseous disk is similar to that of 
Stinson ct al. (2006). Initially, the disk has a simple ex- 
ponential surface density profile. The radial scale length 
of the gas disk, i? ga s, is twice as large as the stellar disk, 
motivated by the observation of Broeils & Rhee (1997). 
We truncate the gas disk at R = 1.5 R gas . Note that the 
truncation radius, i?tnmc = 10-5 kpc, is close to the edge 
of the molecular disk in the Milky Way galaxy (~ 11 kpc; 
Nakanishi & Sofue 2006). The initial vertical distribution 
follows a Gaussian distribution with a scale height of the 
distribution equals to z* (see Table 1). The total gas mass 
is 3.5 x 10 9 Mq. Since we initially place 10 6-7 particles in 
each run, each SPH particle has a mass of 350 — 3500 M e 
(Results of convergence tests are shown in §4.4). The gas 
disk initially rotates with a circular velocity of the model 
galaxy. The gravitational softening length is set to be 
10 pc for the gas particles. The initial gas metallicity is 
set to be the solar value. 

We evolve the disk without radiative cooling for the first 
50 Myr, in order to have a relaxed particle distribution. 
We have checked that our choice of the time duration of 
this initial relaxation phase does not affect the main re- 
sults. 

3.2. Star formation and supernova feedback models 

We adopt a commonly used condition for star forma- 
tion: (1) n H > n th , (2) T < T th , and (3) V • v < 0, for 
a star formation site. We parameterize the star forma- 
tion model by two parameters: n t h and Tth- Here we 
consider two simple models for star formation. One is 
n th = 100 cm" 3 and T th = 5000 K. This density cor- 
responds to mean densities of GMCs, while this tem- 
perature is much higher than the typical temperature of 
GMCs (T < 100 K). Nonetheless, we find that more than 
ninety-percent of stars in mass are formed from the gas 
below T = 100 K. The disk structures are thus insensitive 
to the choice of T t h- We call this model the "high-n t h 
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is the mass of the gas particle, and idyn,i 
t i, respectively. If we use m» iSpawn = m 



gas,u 



1/ - v /47rG/9 g as,i, respectively. If we use 
masses of the gas particles around star forming regions 
become heavier by receiving mass from evolved stars and 
we loose mass resolution. On the other hand, too small 



value of m* 



is not favored from a dynamical point of 



view. We thus fix m* 



to one-third of the original gas 



particle mass as in Okamoto et al. (2003, 2005). When 
the mass of a gas particle becomes smaller than m* jSpawn , 
we convert the gas particle into a collisionless particle. We 
consider each stellar particle as a single stellar population 
(SSP) having its own age and metallicity. We assume the 
Salpeter IMF (Salpeter 1955) whose lower and upper mass 
limits are 0.1 Mq and 100 M , respectively. 

We implement SN feedback in a probabilistic manner 
as in Okamoto ct al. (2008). We assume that stars more 
massive than 8 M Q explode as Type II SNe and each Type 
II SN outputs 10 51 ergs of thermal energy into the ISM 
around the SN. In this paper, we only consider the effect 
of Type II SNe as feedback from stellar particles, since 
the time integration of each run is done only 0.3 — 1 Gyr 
and the lifetime of Type la SNe progenitor is > Gyr. The 
number of SNe in each SSP is approximated by a single 
event. The probability of a SSP i having such event of SN 
explosion during a time interval dt is given by 



PSNII,i 



It 



*SSP,i 
*SSP,i 



\-di 



r SNU (t')dt' 



(7) 



where issp.i is the age of the SSP, rgmi is the SN II rate 
for the SSP, and t$ is the lifetime of a 8 Mq star. SN 
energy is smoothly distributed over the surrounding 32 
SPH particles. We here use the SPH kernel as a weight- 
ing function for the energy deposition. The specific SN 
rate is ~ 0.0072 SN/M© and the typical stellar mass of 
each stellar particle is ~ 1000 M© in our simulation with 
10 6 particles. Thus each SN event in our simulation cor- 
responds approximately to an association of ~ 7 SNe. 

Table 2 shows the model parameters for our runs. We 
perform a total of seven runs, and label them 'A', 'B', 'C, 
'D', and 'E' with additional two runs with extra suffixes, 
such as 3 and 10. Labels indicate the adopted ranges of 
radiative cooling function and star formation model. Run 
A represents a standard model, which is often used in 
cosmological simulations of galaxy formation. This model 
employs 10 4 K as the minimum temperature, T cut . Then 
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Table 1. Parameters of the model galaxy 



DM halo 


Stellar disk 




Gas disk 


M v i T a CNFW b 




A/gas' 


R e R± h 7 ' T- -J 

Jt gas -'Hrunc -^gas lmt 


IO^Mq 12 


4.0 x 10 iu M Q 3.5 kpc 400 pc 


3.5 x 10 y M o 


7 kpc 10.5 kpc 400 pc 10 4 K 



a Virial mass of halo (Mq). b Concentration parameter. c Mass of stellar disk (M©). d Scale length of stellar disk 
(kpc). c Scale height of stellar disk (kpc). Mass of gas disk (Mq). g Scale length of gas disk (kpc). h Truncation 
radius of gas disk (kpc). 'Initial scale height of gas disk (pc). ^Initial temperature of gas disk (K). 



this model does not have the cold phase gas (T < 1000 K). 
In contrast, Run B adopts T cu t = 10 K and has the cold 
phase gas. Both Runs C and D adopt a high-density and 
a low-tempcrature thresholds, whereas these models have 
different star formation efficiencies. Runs C and D em- 
ploy T cut = 10 K. The low star formation efficiency in 
Run C is motivated by the slow star formation model of 
Zuckerman & Evans (1974) and Krumholz & Tan (2007), 
whereas the high star formation efficiency in Run D is mo- 
tivated by the observations of star clusters (C* ~ 0.1 — 0.3) 
reported by Lada & Lada (2003). Such high efficiency 
is also adopted in recent simulations by Tasker & Bryan 
(2006; 2007) (C* = 0.5). Runs C 3 and Ci differ from 
Run C in the mass resolution. Run E excludes the star 
formation process. We evolve Runs A, B, C, D, and E for 
1 Gyr, whereas we terminate Runs C3 and C10 at 0.3 Gyr 
because of limitations of computational resources. 

4. Results 

We present the three-dimensional structure of the ISM 
and the distribution of young stars in §4.1. The SFHs and 
the relation between surface gas density and surface SFR 
density (£ gas — Esfr) are shown in §4.2. Phase (p — T) 
diagrams and density probability distribution functions 
(PDFs) are shown in §4.3. The results of the convergence 
tests arc reported in §4.4. 

4-1- Features of disks 

4-1.1. Structure of gas disks 

Figure 2 shows density and temperature snapshots of 
Runs A, B, C, and D at t = 0.3 Gyr. Different values of 
T cut lead different density and temperature structures; the 
models with T cut = 10 K show gas disks with complex and 
inhomogeneous structures (Runs B, C, and D), whereas 
the model with T cu t = 10000 K show a much smooth gas 
disk (Run A). Different star formation criteria yield fur- 
ther different density and temperature structures. The 
high-nth models have more complex and inhomogeneous 
structures in the ISM (Runs C and D ) than that in the 
low-nth model with T cut = 10 K (Run B). 

Run A has a smooth density structure because the 
higher effective pressure (T cut = 10000 K) stabilizes the 
gas disk. In contrast, Run B has cold clumps because 
of the lower temperature threshold for star formation, 
T cu t = 10 K. Similar influences of T cut on the ISM struc- 
ture is also reported in Saitoh et al. (2006). In Runs A 
and B, the ISM forms stars at local density peaks before 
they develop clumpy structures such as dense filaments, 



which are found in high-nth models. 

Runs C and D have highly inhomogeneous structures 
compared with Runs A and B. In these runs, numerous 
cold and dense (T < 1000 K and nn > 10 cm -3 ) gas clumps 
are formed because of the high-n t h- We can see many 
filaments and clumps of dense gases as well as 'holes' of 
diffuse gases. The dense filaments and clumps have 3 — 4 
orders of magnitudes larger densities than the holes. Such 
structures are rapidly formed after turning on radiative 
cooling and they are retained throughout the evolution. 
Close comparison of Runs C and D reveals that the total 
mass in dense gas clumps is larger in Run C than in Run 
D. The volume which is filled with dense gas in Run C 
also appears larger than in Run D. This is partly because 
the output time for Run D is just after a star formation 
peak (see figure 6). 

From edge-on views of density snapshots in figure 2, we 
find that the gas disks with n-n > 1 cm -3 have a thickness 
of ~ 200 pc for all the models. In Run A, the disk has a 
smooth vertical structure. In Runs B, C, and D, the disks 
appear rather clumpy, similarly to the results of three- 
dimensional Eulerian simulations for galactic disk (e.g., 
de Avillez 2000a; de Avillez 2000b; Tasker & Bryan 2006; 
Wada & Norman 2007). The dense clumps have a low 
temperature (T < 1000 K) in these runs. 

The temperature distributions are shown in the bottom 
row of figure 2. Again, Run A has a smooth and high 
temperature distribution, while Runs B, C, and D have 
complex structures of cold gas. It is intriguing that cold 
components with T < 1000 K in Runs B, C, and D are not 
always associated with dense gas layers (see close up views 
of density and temperature maps from edge-on). Dense 
gas layers are heated up by hydrodynamic shocks caused 
by gravity and SNe. As a result, the temperature at a 
given density can be different from the equilibrium value. 
This result is similar to that of previous ISM simulations 
of galactic disks and circumnuclear disks (e.g., Wada & 
Norman 1999; Wada & Norman 2001; Wada & Tomisaka 
2005; Tasker & Bryan 2006). 

We define a characteristic scale height for each phase 
as the height that contains half of the mass of the phase 
at given radius. We consider two phases, a cold phase 
(10 K < T < 100 K) and a warm phase (100 K < T < 
10000 K) and call their characteristic scale heights, z co id 
and Zwarm, respectively. Figure 3 shows z co id (Runs B, C, 
and D) and z warm (Runs A, B, C, and D) as a function of R 
at t = 0.3 Gyr. Immediately, we find that typical vertical 
distribution is clearly separated by gas temperature and 
the distribution is not affected by criteria of star formation 
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Fig. 2. Density (top row) and temperature (bottom row) snapshots of Runs A, B, C, and D (from left to right), at i = 0.3 Gyr. 
Top and middle panels show edge-on views (a thin slice of y = 0) while the bottom panel shows a face-on view (a thin slice of z = 0). 
The middle panels show a region of —15 < x < 15 kpc and —2.5 kpc < z < 2.5 kpc. The top panels show a region of < x < 5 kpc 
and —0.42 kpc < z < 0.42 kpc. The plot range of the bottom panels is —15 < x < 15 kpc and —15 < y < 15 kpc. The insets in the 
bottom panels give the close up view of the inner disk for the first quadrant (0 < x < 5 kpc and < y < 5 kpc). 
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Table 2. Parameters of runs 



model 


N a 


^SPH b 


e 


c 


rp d 

J cut 


rath 6 




TV 


c* 


Run A 


10 b 


3500 M Q 


10 


pc 


10000 K 


0.1 cnT 




15000 K 


0.033 


Run B 


10 6 


3500 M Q 


10 


pc 


10 K 


0.1 cm" 


-3 


15000 K 


0.033 


Run C 


10 6 


3500 M 


10 


pc 


10 K 


100 cm" 


-3 


5000 K 


0.033 


Run C 3 


3x 10 6 


1170 M Q 


10 


pc 


10 K 


100 cm" 


-3 


5000 K 


0.033 


Run Cio 


10 7 


350 Mq 


10 


pc 


10 K 


100 cm" 


-3 


5000 K 


0.033 


Run D 


10 6 


3500 M Q 


10 


pc 


10 K 


100 cm" 


-3 


5000 K 


0.5 


Run E 


10 6 


3500 M Q 


10 


pc 


10 K 


N/A 




N/A 


N/A 



a The initial number of SPH particles. b Mass of individual SPH particles (Mq). c Gravitational softening length (pc). 
d Cut off temperature of cooling function (K). °Threshold density of star formation (cm -3 ). f Thrcshold temperature 

of star formation (K). g Star formation efficiency. 
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Fig. 3. Half mass heights of gas disks with T < 100 K and 
100 K < T < 10000 K at t = 0.3 Gyr. The thick lines represent 
2-coldj while the thin curves represent ^ warm . Solid (black), 
dashed (red), dot-dash (green), and dash-dot-dot-dot (blue) 
curves indicate Runs A, B, C, and D, respectively. The scale 
of cold gas disk (z co i<) ) for Run A is not plotted on this figure. 

in our simulations. In Runs B, C, and D, z co i d is 20 — 50 pc, 
whereas x warm is 60 — 120 pc. Run D has two local peaks in 
the profiles: one is shown in the curve of z warm at i?~ 5 kpc 
and another one is shown in that of z co i d at R > 8 kpc. 
These two peaks are induced by the localized SNe and are 
temporary structure. In Run A, z warm is identical to the 
others, except at the large radii (R > 7 kpc). Hence the 
vertical structures of the ISM are basically determined by 
the ranges of the cooling function. Both z co id and z waim 
in all runs gradually increase with increasing R. In Run 
A, the curve of z W arm becomes four times thicker at the 
edge than that at center. In Runs B, C, and D, the curves 
of Zcoid and z warm become two times thicker at the edge 
than those at centers. 
4-. 1.2. Structure of stellar disks 

Figure 4 shows the face-on and edge-on views of the 
distribution of star particles for the four runs. We color 
the figures based on the age of the ten-percent youngest 
star particles in each grid as a representative value. The 
face-on views clearly show that the distributions of stars 
in Runs C and D have clumpy structures while Runs A 
and B have smooth structures. The localization of young 
stars in high-nth runs is easily understood by the complex 



and inhomogeneous structure of the gas disk. Close com- 
parison between figures 2 and 4 reveals that young stars 
are not always associated with dense gas regions. This 
is because SN feedback blows out remaining gas around 
young stars. Consequently, no clear spatial correlation is 
found between the distributions of young stars and that 
of dense gas clumps at any given time. 

The edge-on views in figure 4 suggest that the verti- 
cal distribution of star particles is strongly affected by 
the threshold density for star formation. This result indi- 
cates that nth is an essential parameter that determines 
the thickness of stellar disk. We argue that the threshold 
density needs to be determined by the physical value of 
the star forming regions, such as typical density of molec- 
ular clouds. 

We further study the vertical structures of the stellar 
disks. We define characteristic scale heights as the heights 
that contains 50%, z*,5o, and 90%, z*,9o, of the stellar 
mass at a given radius. Figure 3 shows z gas ,50 an d -Zgas,90 
as a function of R at t = 0.3 Gyr. This indicates clearly 
that the threshold density for the star formation strongly 
affects the vertical structure of stellar disks. The low-n t h 
models (Runs A and B) have thick, extended stellar disks 
(30— 110 pc for z* j50 and 100 — 300 pc for 2*,go), whereas 
the high-nth models (Runs C and D) have very thin stellar 
disks (10 — 30 pc for z*. 50 an d 30 — 60 pc for z* t 9o). We also 
find that heights of all stellar disks increase, even weakly, 
with R. The threshold density has a critical role for stellar 
disk formation. In contrast, we note that the value of 
C* does not significantly affect the vertical structures of 
stellar disks. 

4- 1-3. Disk scale heights: comparison of simulations and 
observations 

In this subsection, we compare the vertical scale heights 
of three components in our models with observations, 
namely the vertical distributions of HI and H2 gas disks, 
and young star forming regions (galactic young open clus- 
ters). To simplify the comparison, we utilize half-mass 
heights both simulations and observations. 

First, we compare z warm with the half-mass height of 
the galactic HI gas disk, since the typical temperature of 
HI gas is 100 K < T < 10000 K (Myers 1978; Spitzer 1978). 
Nakanishi & Sofuc (2003) reconstructed the three dimen- 
sional structure of HI gas of the Milky Way galaxy, by 



8 Saitoh et al. [Vol. , 

10 7 1° 6 1° 9 Age [year] 



x \kpc] y. \kpc] x fkpcl y. [kpcl 

1 2 J 4 G 1 2 T 4 G 1 2 T 4 1 2 J 4 5 




j I I I I I I I I I I I I I ill | | | | | | L *i ■ I | | | | ill | | | | | | | | | | | | | l±j I I I I I I I I I I I I I L 

-10 10 -10 10 -10 10 -10 10 
k [kpc] Run A x [kpc] Run B x [kpc] Run G * [kpc] Run D 



Fig. 4. Projected star particle distributions in face-on and edge-on views of Run A, B, C, and D. Star particles are assigned on 
uniform grids. Each grid size is 234 pc = (30 kpc/128) for for middle and bottom panels, whereas that is 39 pc for top panels and 
the insets in bottom panels. The color level calculated based on the age of the ten-percent youngest star particles in each grid. 
Arrangement of panels is the same as figure 2. 
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Fig. 5. Half mass and ninety-percent mass heights of 
stellar disks at t = 0.3 Gyr. The thick lines repre- 
sent z*,s0i while the thin curves represent z*,go. Solid 
(black), dashed (red), dot-dash (green), and dash-dot-dot— 
dot (blue) curves indicate Runs A, B, C, and D, respectively. 



compiling three HI survey data: the Leiden/Dwingelloo 
survey (Hartmann & Burton 1997), Parkes survey (Kerr 
et al. 1986), and NRAO survey (Burton & Liszt 1983). 
They obtained the vertical scale height, which is de- 
fined as the full width at half maximum (FWHM) as 
the vertical scale height, as a function of R (see figure 
4 in Nakanishi & Sofue 2003). To multiple 0.625 for the 
FWHM, we obtain the half-mass scale height, since they 
assume the vertical distribution of HI gas is proportional 
to the square of a function of hyperbolic secant. The ob- 
servationally suggested half-mass scale height (60— 180 pc 
at R = — 10 kpc) is almost identical to z warm for the four 
runs (Runs A, B, C, and D). 



Second, we compare z co id and the half-mass scale height 
of the galactic H2 gas disk. By using the compilation 
data of 12 CO( J =1-0), which is provided by Dame et al. 
(2001), Nakanishi & Sofue (2006) found that the scale 
heights of the H 2 disk (the FWHM of the H 2 gas disk) in 
the Milky Way galaxy is 48 - 160 pc at R = - 11 kpc. 
Again, we convert the FWHM to the half mass scale. The 
half mass scale of the observed H 2 gas disk is ~ 30 — 100 pc 
at R = — 11 kpc. The half-mass scale heights of simu- 
lations, z co id, are slightly thiner than that in the Milky 
Way galaxy. However the difference is at most a factor 
of two. Then we consider that z C oid is in good agreement 
with that in the Milky Way galaxy. 

Finally, we compare the vertical scales of star particles 
with the observed one through mass fractions. The verti- 
cal distribution of galactic young open clusters is a good 
tracer of the recent star forming regions in the Milky Way 
galaxy. As shown in Janes & Phelps (1994), the verti- 
cal distribution of galactic young clusters of which the 
ages are shorter than that of Hyades (~ 800 Myr) is fit- 
ted by an exponential function with a 55-pc scale- height. 
The half-mass height of the exponential profile is ~ 38 pc 
(0.69 x 55 pc). Star particles in high-n t h models (Runs C 
and D) have z*,5o ^10 — 30 pc, whereas stars in the low- 
nth models (Runs A and B) are more broadly distributed 
in the vertical direction, z».5o ~ 60 — 100 pc. Only the 
high-nth models are consistent with the observation. It 
should be noted that even in Run B which includes the 
cooling under 10 4 K, the scale height is too large. We 
can conclude that n t h is a key parameter to determine the 
thickness of stellar disks. 




Time [Gyr] 

Fig. 6. Star formation histories in the simulations. Solid, 
dashed, dot-dash, and dash-dot-dot-dot lines indicate Runs 
A, B, C, and D, respectively. 

4-2. SFHs and S gas — Xsfr relations 

Figure 6 compares the global SFHs in four runs (Runs 
A, B, C, and D). In all runs, the SFHs are character- 
ized by an initial rapid increase, followed by a gradual 
decrease. The rapid increase continues for only ~ 10 7 yr 
because of SN feedback effects. The evolutions after the 
first peak are described approximately by an exponential 
decay with short-period spikes. We find that the SFRs 
asymptotically decrease to ~ 1 — 2 M Q yr ~ 1 , which is close 
to the observational value of SFR in nearby spiral galax- 
ies (~ 1 M yr -1 : James et al. 2004) and Galactic SFR 
(~ 3 M Q yr" 1 : McKee & Williams 1997). 

It is worth pointing out that the difference in the global 
SFHs of Runs C and D is rather small, despite a factor of 
15 (= 0.5/0.033) difference in C*. Tasker & Bryan (2006) 
performed test simulations by changing C* by a factor of 
ten and found that it has little effect to the global SFHs. 
In other words, global SFHs are not directly proportional 
to C* when a high-density threshold is adopted. They 
claimed that this is because the dynamical time in star 
forming regions is sufficiently shorter than that in galaxies. 
As will be discussed in §4.3.2, we find that the contraction 
timescale of the gas is about 5 times longer than the local 
dynamical time and this timescale does not depend on 
the value of C*. Hence the global SFH is not directly 
proportional to C* . 

Figure 7 presents S gas — Esfr relations for Runs A, B, 
C, and D. All of our runs show similar relations, both 
in the slope and in the normalization, to the observa- 
tional values (-1.4; solid lines in the figure), although 
the range covered by our simulations is somewhat lim- 
ited. Slightly steeper slopes in Runs with a high critical 
density (i.e., Runs C and D) than Run A is consistent with 
what Wada & Norman (2007) found in their theoretical 
model of global star formation in the ISM with log-normal 
density PDF. We have shown that the different threshold 
models have very different ISM structures. The distri- 



Fig. 7. Surface gas density and the surface star formation 
rate for Runs A, B, C, and D in three different epochs 
t = 0.3, 0.5, and 1.0 Gyr. The surface SFRs are computed us- 
ing the surface densities of young star particles of which ages 
are shorter than the typical age of massive stars; 4.5 X 10 7 yr. 
The outer edges of surface densities correspond with the dis- 
tance of the most distant young star particle from the galac- 
tic center and the typical edges of the star forming regions 
are R ~ 10 kpc. Four cylindrically averaged values with a 
constant radial interval are obtained from each run and each 
epoch. Stars, crosses, circles, and squares represent the se- 
quences of Runs A, B, C, and D, respectively. The solid line 
is a best fit from observations (Kennicutt 1998). 

butions of newly formed star particles are also different, 
reflecting the structures of the ISM. Nevertheless, as we 
have shown, the observed E gas — Ssfr relation is well re- 
produced in all our simulations. 

4-3. Phase structure of the ISM 

In this section, we study the phase structures of the ISM 
in more detail. The differences between the low-nth and 
high-nth models in phase diagram are studied in §4.3.1. 
The distribution functions of gas mass as a function of 
density and temperature are also discussed. Evolution 
of individual gas particles is analyzed in §4.3.2. We find 
from the analysis that the evolution timescale of the ISM 
(1 cm -3 < ?ih < 100 cm~ 3 ) is typically ~ 5 idyn(nH")- We 
show the PDFs of Runs C, D and E in §4.3.3. The high- 
density parts of PDFs in Runs C and D are well fitted by 
a log-normal function, whereas that in Run E is a single 
power-law like form. In our simulations, a form of PDF 
changes by the effects of star formation and SN feedback. 
4-3.1. Phase diagram 

In figure 8, we show global phase (p — T) diagrams of 
Runs A, B, and C at t = 0.3 Gyr. The phase diagram 
in Run D resembles that in Run C and thus we exclude 
it here. Differences in phase distribution between Runs 
C and D are discussed in §4.3.3. Phase structures in dif- 
ferent runs are very different. Phase structures strongly 
depend on n t h and T cut . As a consequence of the assump- 
tion (T cut = 10 4 K), Run A consists of a warm (T~ 10 4 K) 
and a hot (T > 10 5 K) phase. There is no cold phase be- 
cause of the cut-off temperature of the cooling function, 
T C ut, whereas the hot phase is formed by the SN feed- 
back. In Run A, there is a density cut off at ~ 10 cm -3 
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to 



Pin /(A Pm /At) 



(8) 



where p m is initial values ol median densities and Ap m 
is density changes of p m within At, respectively. We 
then find that the evolution time is t cvo ~ 5 idyn^H): 
^cvo("-h ~ 100 cm -3 ) ~ 20 Myr and t cvo (nn ~ 1 cm -3 ) ~ 
130 Myr. The analysis on Run D shows almost the same 
results (compare the thin (Run D) and solid (Run C) 
lines in figure 10). Interestingly, the density dependence 
is proportional to the local dynamical time: t cvo oc p^ 1 ^ 2 . 
Comparison between Runs C and D indicates that the 
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because star formation is rapidly consuming high-density 
gas and the gas pressure with T = 10 4 K prevents further 
collapse. The overall feature is similar to that in Stinson 
et al. (2006) and the assumption (T cut = 10 4 K) is reason- 
able for coarse resolution (mgpjj ~ 10 5-6 M©) runs (see 
figure 1). Runs B and C show a clear multiphase structure 
with a cold gas. The high-density gas in Run B becomes 
denser than that in Run A due to the low effective pres- ^ 
sure. The high-density tail extends to ~ 100 cm~ 3 in this > — > 
case. The high-nth model has more dense and cold gas, | — 
since the gas consumption due to star formation is oc- 
curred at ~ 100 cm -3 . The multiphase structure in Run 
C is similar to those obtained in previous high-resolution 
grid-base simulations (e.g., Wada & Norman 1999; Wada 
2001; Tasker & Bryan 2006). 

The mass fractions as a functions of density F(tih) and 
temperature F(T) are shown in top and right histograms 
in each panel of figure 8. In Run A, the peak of F(nn) 
is n H ~ 1 cm" 3 and that of F(T) is T ~ 10 4 K. The 
temperature peak is shifted to ~ 100 K in Run B. The 
dominant component in temperature in Run C is also the 
cold phase gas (T < 1000 K) as is also found in other high- 
resolution simulations of the ISM, i.e., ~80 — 90% of mass 
in the ISM is in the cold phase (i.e., Rosen & Bregman 
1995; Wada 2001; Tasker & Bryan 2006). The peak of 
F(nn) is n H ~ 1 cm" 3 and that of F(T) is T ~ 100 K. 
In conclusion, the multiphase ISM in the high-n t h models 
has a large mass of reservoir around 1 cm -3 for the star 
forming region (nn > 100 cm -3 ). Thus the evolution of 
the reservoir should play a key role for the star formation 
in the multiphase ISM. We further study the evolution of 
the gas from the reservoir to star forming regions in the 
next subsection. 

4-3.2. Detailed evolution of fluid elements in the multi- 
phase ISM 

Figure 9 shows the evolution of gas in the phase (p — T) 
diagram within At = 1 Myr at t = 0.3 Gyr in Runs A and 
C. There is a main stream toward higher densities between 
T cq and 10° 5 x T cq in Run C, while we can not find any 
flow in that direction in Run A. We compute the mass flux 
on the phase diagram around the main stream, between 
T eq and 10 x T eq , toward high density with binning from 
log(rin) = to log(nn) = 2 every 1/3 dex in density. The 
median values are adopted for the indicator of the density 
evolution in each bin. Figure 10 shows the e-folding (evo- 
lution) times in Runs C and D. We define the evolution 
time that the density changes e times by the evolution 
time of medians within At interval on the phase diagram: 
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Fig. 8. p — T diagram for Runs A, B, and C (from the 
top panel to the bottom panel). Whole plot regions are 
subdivide into 128 X 128 grids and each grid is colored by 
a mass fraction. Top and right histograms in each panel 
show mass fractions as functions of density and temperature. 
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Fig. 9. Evolution of SPH particles in each grid on phase 
diagram within At = 1 Myr at t = 0.3 Gyr in Runs A (top 
panel) and C (bottom panel). Arrows indicate evolution 
of sampled grids on the phase diagram in T < 10 5 K for 
Run A and in T cq < T < 10 X T oq for Run C. The heads 
of arrows indicate the median values of density and tem- 
perature after At of evolution. Black curves indicate the 
equilibrium temperature of the adopted cooling and heat- 
ing with the solar abundance, T eq , and its families that 
10 ' 5 T eq and 10 T eq , respectively. Background gray scales 
indicate mass weighted distributions of gas on phase diagram. 

evolution timescale in the multiphase ISM is independent 
of the star formation efficiency, C* . 

Figure 11 shows density evolutions of several selected 
SPH particles in Run C. There are five loci of randomly 
selected SPH particles from 190 Myr to 280 Myr. All of 
the particles have the density changes of ~ 5 orders of 
magnitudes. The evolutions are not monotonic but com- 
plex with many compressions and expansions. The con- 
siderable drive mechanisms of the compressions are grav- 
itational collapses and hydrodynamical shocks, whereas 
those of the expansions are SNe and shear extensions. The 
global evolution of the ISM is expressed by the superpo- 
sitions of above mechanisms. Hence the mass flow from 



Fig. 10. The e-folding times of density evolution for Run 
C (the thick solid line with circles) and Run D (the 
thin solid line with boxes). See the text for the defi- 
nition of the evolution timescale. Dotted line represents 
the local dynamical time as a function of the density. 
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Fig. 11. Evolution tracks of densities for randomly se- 
lected SPH particles from 190 Myr to 280 Myr in Run C. 
Five colors denote five different particles. Corresponding lo- 
cal dynamical time for density is shown in the right side. 



~ 1 cm 3 to ~ 100 cm 3 is decelerated compared with 
that expected by only the gravitational collapse. 
4.3.3. The density PDF 

Shapes of PDF of the ISM give us an idea to con- 
nect the global structure of the ISM and star forming re- 
gions. There have been a number of studies on this issue 
(e.g., Vazquez- Scmadcni 1994; Scalo et al. 1998; Vazqucz- 
Semadeni et al. 2000; Wada 2001; de Avillez & Mac 
Low 2002; Kravtsov 2003; de Avillez & Breitschwerdt 
2004; Slyz et al. 2005; Wada & Norman 2007; Tasker 
& Bryan 2007; Robertson & Kravtsov 2007). When the 
fluid evolution is dominated by a random density change 
without any scale dependence, it is expected that the 
volume weighted PDF becomes 'log-normal' distribution 
(Vazquez-Semadeni 1994). In the multiphase ISM, there 
are many physical processes such as radiative cooling, 
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FUV heating, star formation, SN feedback, hydrodynam- 
ical shock, and self-gravity. It is unclear that whether the 
combination of these processes is either scale-invariant or 
not. The expected form of a resultant PDF is also unclear. 

Figure 12 shows the evolutions of PDFs in our SPH 
simulations. Here, we plot only PDFs of Runs C, D, and 
E. The evolutions of PDFs in these three runs are almost 
the same at the first phase, although the quasi-static fi- 
nal shapes depend on models. At initial state, the PDFs 
have uniform density distributions with a peak around 
riH ~ 0.3 cm -3 . The PDFs are smoothed out within 
first ~ 0.2 Gyr. Our PDFs in the three runs then be- 
come quasi-static as reported by Wada & Norman (2001), 
Kravtsov (2003), and Wada & Norman (2007). High den- 
sity tails (riH > 1 cm~ 3 ) in Runs C and D are in steady 
state and have log- normal forms (e.g., Vazquez-Semadeni 
1994; Vazquez-Semadeni et al. 2000; Wada 2001; Kravtsov 
2003; Wada & Norman 2007; Tasker & Bryan 2007). The 
high-density tail in Run E, the model does not include 
cither star formation or SN feedback, is asymptotic to a 
power-law form (green and blue lines). The difference in 
high-density region comes from the fact that the gas in the 
high-density region converts into stars and SN feedback 
blows out the surrounding dense gas in Runs C and D. 
Slyz et al. (2005) showed that a high-density tail of a PDF 
has a power-law like form in a self-gravity-dominatcd sys- 
tem, although several authors argued that the log-normal 
PDF is the robust structures of multiphase ISM, regard- 
less of the input physics (Wada & Norman 2001; Kravtsov 
2003) . Further investigation into details is required for the 
response of the input physics (e.g., UV radiation, Susa & 
Wada in preparation) in the form of PDF. 

4-4- Convergence tests 

Figure 13 shows density maps in the simulations us- 
ing N = 10 6 ,3 x 10 6 , and 10 7 (Runs C, C 3 , and Ci ) at 
t = 0.3 Gyr. At the time, the numbers of particles in 
all of runs increase ~ 30 % from the initial states due to 
star formation. The density distributions indicate that 
differences among the three runs with different mass res- 
olutions are very little. These gas disks have almost the 
same sizes of filaments and voids. The complcxness of 
the gas disks is almost the same degree. The statistical 
structures of the multiphase ISM, PDFs, in the three runs 
are also similar to one another (see the lower-right panel 
in figure 13). Differences are found in the high-density 
tails from residuals. However the residuals are only factor 
of three at most. Figure 14 shows vertical thicknesses of 
stellar disks. These thickness are almost the same for all 
runs. Thus we consider that our results are roughly con- 
verged in above mass resolutions, for our selected values 
for n t h and T cut . 

5. Summary and Discussion 

5.1. Importance of the density threshold for star forma- 
tion 

In many studies, numerical models of star formation 
in galaxy formation have been calibrated to be consistent 
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Fig. 12. Volume weighted probability distribution functions 
(PDFs) of Runs C, D, and E (from top to bottom) within 
-R < 10 kpc. We calculate a volume fraction of an SPH par- 
ticle by Vi = vail pi. The density interval is 0.2 dex. The 
PDFs are normalized to unity. Runs C and D show six differ- 
ent epochs (t = 0, 0.1, 0.2, 0.3, 0.5, and 1.0 Gyr), while 
Run E shows first four epochs. Black, red, green, blue, 
sky blue, and magenta indicate PDFs of different epochs 
t = 0, 0.1, 0.2, 0.3, 0.5, and 1.0 Gyr, respectively. The ver- 
tical dotted lines are corresponding with the star formation 
threshold density, n t h = 100 cm -3 . The high-density region 
(njj > 100 cm -3 ) in Run E is hatched, since the region has 
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Fig. 13. Sets of panels at tops and left-bottom are gas density maps for three different mass resolu- 
tions (JV = 10 6 , 3 X 10 6 , and 10 7 , respectively). The right-bottom panel shows PDFs of three runs at 
t = 0.3 Gyr. Thin dashed, thick solid, and thin solid lines indicate simulations using N = 10 6 , 3 X 10 6 , and 
10 7 , respectively. PDF residuals of Runs C3 and C10 from C is shown in the upper portion of the panel. 
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with the observational £ gas — Esfr relation. For example, 
Stinson et al. (2006) successfully reproduced the relation 
with the threshold density of n t h = 0.1 cm -3 . In order 
to reproduce the relation within simulated disk galaxies, 
Springel & Hernquist (2003) and Schaye & Dalla Vecchia 
(2007) explicitly involved the Schmidt-Kcnnicutt relations 
in star formation models in the ISM. Kravtsov (2003) 
showed that a high-density threshold (n t h = 50 cm -3 ) with 
a constant star formation time (= 4 Gyr) reproduce the 
S gas - E SF R relation. Tasker & Bryan (2006, 2007) re- 
ported the comparison of the low-nth and high-n t h star 
formation model. Their simulations revealed that both 
models are able to reproduce the observed £ gas — Ssfr 
relation. We examined two models, a low-nth model 
(nth = 0.1 cm -3 : Runs A and B) and a high-n t h model 
(«th = 100 cm -3 : Runs C and D). Our results also show 
that both models can reproduce the E gas — Xsfr relation 
(see figure 7). 

We argue that we should generate "stars" above the 
physical density of real star forming regions such as GMCs 
or molecular cores to investigate the detailed structure 
and evolution of disk galaxies. In this paper, we high- 
lighted on the ISM structure and the distribution of newly 
formed star particles. We found that (1) only the high- 
nth model reproduces the complex, inhomogcncous, and 
multiphase ISM structures (see figures 2 and 4), where the 
cold gas dominates in mass (see figure 8). These natures 
are well comparable with other studies of the ISM: the 
geometrically complex and inhomogcncous structures of 
the ISM (e.g., Rosen & Bregman 1995; Wada & Norman 
1999; de Avillez 2000a; de Avillez & Berry 2001; Tasker 
& Bryan 2006; Tasker & Bryan 2007), three phases struc- 
tures of the ISM where the cold mass dominates (e.g., 
McKee & Ostriker 1977; Myers 1978; Rosen & Bregman 
1995; Wada 2001; Tasker & Bryan 2006; Tasker & Bryan 
2007). The log- normal PDF in the ISM is also found in 
the high-nth models (see figure 12 and figure 13) although 
the origin of the log-normal shape appear to be differ- 
ent from that in previous studies (e.g., Vazqucz-Scmadcni 
1994; Vazquez- Semadeni et al. 2000; Wada 2001; Kravtsov 
2003; Wada & Norman 2007; Tasker & Bryan 2007). (2) 
Only the high-nth models can reproduce observationally 



reported scale heights of gas disks and young star forming 
regions as shown in figures 3 and 5. Therefore we empha- 
size that the density threshold for star formation model 
is the key parameter to model realistic three-dimensional 
structures of galaxies, especially gas and stellar disk struc- 
tures. We have to choose the star forming gases as the 
physical one for this purpose. It is necessary to solve en- 
ergy equation for much lower temperature gas than 10 4 K 
to resolve the high-density gas. This requires higher mass 
resolution than those used in previous simulation of galaxy 
formation. 

5. 2. Weak dependence on star formation efficiency 

Runs C and D differ in the values of C*. The results are, 
however, similar in terms of the ISM structure and stellar 
disks (figures 2 and 4), the SFHs (figure 6), the E gas — 
Ssfr relation (figure 7), and the phase distribution of the 
ISM (figure 12). Why do the simulations show similar 
results? As shown in figure 8, a large fraction of the gas 
exist at around ~ 1 cm~ 3 and it behaves as the reservoir of 
the star forming gas. The mass supply timescale from the 
reservoir to the star forming region determines the global 
star formation rate in the model that adopts the high- 
density threshold (n t h = 100 cm -3 ). From figure 10, we 
find that the timescale is ~ 5 idyn(^H) and the timescale is 
independent of the values of C* . Hence C* has only weak 
effects on the global features of the ISM and stellar disks 
formed from the ISM. 

One of the most important advantages of the high-nth 
model is that global SFR is not determined by the local 
quantities but by the global state of the ISM. The exact 
value of C* is unimportant in galactic scale simulations, 
hence we can effectively avoid an uncertainty in the star- 
formation model. The remaining weak dependences of C* 
would vanish when we adopt even higher threshold density 
for star formation together with higher resolution. 

When we tune the value of C* in low-nth models, we 
can obtain the similar global star formation properties 
in disk galaxies regardless of whether we resolve the de- 
tailed structure of the ISM or not. This finding provides 
a support for star formation models in lower resolution 
cosmological simulations, which is similar to Run A in 
this paper. We will investigate whether it is also true for 
starburst galaxies in forthcoming papers. 
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